www.gusucode.com > C++ 图像分割源代码 > C++ 图像分割源代码/gusucode/seg_source/Analysis.cpp
//Download by http://www.NewXing.com //Copyright (c) 2004-2005, Baris Sumengen //All rights reserved. // // CIMPL Matrix Performance Library // //Redistribution and use in source and binary //forms, with or without modification, are //permitted provided that the following //conditions are met: // // * No commercial use is allowed. // This software can only be used // for non-commercial purposes. This // distribution is mainly intended for // academic research and teaching. // * Redistributions of source code must // retain the above copyright notice, this // list of conditions and the following // disclaimer. // * Redistributions of binary form must // mention the above copyright notice, this // list of conditions and the following // disclaimer in a clearly visible part // in associated product manual, // readme, and web site of the redistributed // software. // * Redistributions in binary form must // reproduce the above copyright notice, // this list of conditions and the // following disclaimer in the // documentation and/or other materials // provided with the distribution. // * The name of Baris Sumengen may not be // used to endorse or promote products // derived from this software without // specific prior written permission. // //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT //HOLDERS AND CONTRIBUTORS "AS IS" AND ANY //EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT //NOT LIMITED TO, THE IMPLIED WARRANTIES OF //MERCHANTABILITY AND FITNESS FOR A PARTICULAR //PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE //CONTRIBUTORS BE LIABLE FOR ANY //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, //EXEMPLARY, OR CONSEQUENTIAL DAMAGES //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT //OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, //DATA, OR PROFITS; OR BUSINESS INTERRUPTION) //HOWEVER CAUSED AND ON ANY THEORY OF //LIABILITY, WHETHER IN CONTRACT, STRICT //LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR //OTHERWISE) ARISING IN ANY WAY OUT OF THE USE //OF THIS SOFTWARE, EVEN IF ADVISED OF THE //POSSIBILITY OF SUCH DAMAGE. #include "./Analysis.h" #include "mkl.h" namespace Analysis { void StatusDisplay(long status) { long classError; classError = DftiErrorClass(status, DFTI_ERROR_CLASS); if(! classError) { cerr << "Error Status is not a member of Predefined Error Class\n"; } else { char* errorMessage = DftiErrorMessage(status); cerr << "Error_message = " << errorMessage << endl; } } Vector<ComplexFloat> FFT(Vector<ComplexFloat>& x) { int n = x.Length(); Vector<ComplexFloat> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Vector<ComplexFloat>& FFTI(Vector<ComplexFloat>& x) { int n = x.Length(); //Vector<ComplexFloat> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } Vector<ComplexDouble> FFT(Vector<ComplexDouble>& x) { int n = x.Length(); Vector<ComplexDouble> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Vector<ComplexDouble>& FFTI(Vector<ComplexDouble>& x) { int n = x.Length(); //Vector<ComplexDouble> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } Vector<ComplexFloat> FFT(Vector<float>& x) { int n = x.Length(); Vector<ComplexFloat> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_REAL, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } for(int i=1;i<(n+1)/2;i++) { y[n-i] = conj(y[i]); } return y; } Vector<ComplexDouble> FFT(Vector<double>& x) { int n = x.Length(); Vector<ComplexDouble> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_REAL, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } for(int i=1;i<(n+1)/2;i++) { y[n-i] = conj(y[i]); } return y; } // //////////////////// // Backward transform // //////////////////// Vector<ComplexFloat> IFFT(Vector<ComplexFloat>& x) { int n = x.Length(); Vector<ComplexFloat> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Backward transform status = DftiComputeBackward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Vector<ComplexFloat>& IFFTI(Vector<ComplexFloat>& x) { int n = x.Length(); //Vector<ComplexFloat> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Backward transform status = DftiComputeBackward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } Vector<ComplexDouble> IFFT(Vector<ComplexDouble>& x) { int n = x.Length(); Vector<ComplexDouble> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Backward transform status = DftiComputeBackward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Vector<ComplexDouble>& IFFTI(Vector<ComplexDouble>& x) { int n = x.Length(); //Vector<ComplexDouble> y(n); DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Backward transform status = DftiComputeBackward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } // ///////////////// // 2D Transform // ///////////////// Matrix<ComplexFloat> FFT2(Matrix<ComplexFloat>& x) { int rows = x.Rows(); int cols = x.Columns(); Matrix<ComplexFloat> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Matrix<ComplexFloat>& FFT2I(Matrix<ComplexFloat>& x) { int rows = x.Rows(); int cols = x.Columns(); //Matrix<ComplexFloat> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } Matrix<ComplexDouble> FFT2(Matrix<ComplexDouble>& x) { int rows = x.Rows(); int cols = x.Columns(); Matrix<ComplexDouble> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Matrix<ComplexDouble>& FFT2I(Matrix<ComplexDouble>& x) { int rows = x.Rows(); int cols = x.Columns(); //Matrix<ComplexDouble> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeForward( DescHandle, x.Data() ); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } Matrix<ComplexFloat> IFFT2(Matrix<ComplexFloat>& x) { int rows = x.Rows(); int cols = x.Columns(); Matrix<ComplexFloat> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols)); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeBackward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Matrix<ComplexFloat>& IFFT2I(Matrix<ComplexFloat>& x) { int rows = x.Rows(); int cols = x.Columns(); //Matrix<ComplexFloat> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols)); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeBackward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } Matrix<ComplexDouble> IFFT2(Matrix<ComplexDouble>& x) { int rows = x.Rows(); int cols = x.Columns(); Matrix<ComplexDouble> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols)); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeBackward( DescHandle, x.Data(), y.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return y; } Matrix<ComplexDouble>& IFFT2I(Matrix<ComplexDouble>& x) { int rows = x.Rows(); int cols = x.Columns(); //Matrix<ComplexDouble> y(rows,cols); int dims[2] = {cols,rows}; int strides[3] = {0,rows,1}; DFTI_DESCRIPTOR_HANDLE DescHandle = 0; long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Strides status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Not Inplace //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //if(status != DFTI_NO_ERROR) //{ // StatusDisplay(status); // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Problem with FFT!"); //} // Backward scale status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols)); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Commit Dfti descriptor status = DftiCommitDescriptor(DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } // Compute Forward transform status = DftiComputeBackward( DescHandle, x.Data()); if(status != DFTI_NO_ERROR) { StatusDisplay(status); cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Problem with FFT!"); } status = DftiFreeDescriptor(&DescHandle); if(status != DFTI_NO_ERROR) { StatusDisplay(status); Utility::Warning("Problem when trying to free the FFT descriptor handle!"); } return x; } //========================================================================= // Convolution //========================================================================= Vector<float> Conv(Vector<float>& input, Vector<float>& filter, FilterArea fa, Border b, bool PowerOfTwo) { return Real( Conv(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) ); } Vector<double> Conv(Vector<double>& input, Vector<double>& filter, FilterArea fa, Border b, bool PowerOfTwo) { return Real( Conv(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) ); } Vector<ComplexFloat> Conv(Vector<ComplexFloat>& input, Vector<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo) { // Check if the filter size is smaller than the input size if(filter.Length() > input.Length()) { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Filter length should not be larger than the input signal!"); } int fftlength = input.Length()+filter.Length()-1; if(PowerOfTwo) { fftlength = NextPow2(fftlength); } Vector<ComplexFloat> temp1; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp1 = Vector<ComplexFloat>(fftlength); for(int i=0; i<input.Length(); i++) { temp1[i] = input[i]; } } else if(b == Circular) { temp1 = input; } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } if(b == Symmetric) { int wid1 = filter.Length()/2; int wid2 = filter.Length() - wid1 - 1; if(wid1 != 0) { Vector<ComplexFloat> c1 = Reverse( input.Slice(0, wid1-1) ); temp1.ReadFromVector(c1, temp1.Length()-wid1); } if(wid2 != 0) { Vector<ComplexFloat> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) ); temp1.ReadFromVector(c2, input.Length()); } } else if(b == Replicate) { int wid1 = filter.Length()/2; int wid2 = filter.Length() - wid1 - 1; if(wid1 != 0) { Vector<ComplexFloat> c1(wid1, input.ElemNC(0)); temp1.ReadFromVector(c1, temp1.Length()-wid1); } if(wid2 != 0) { Vector<ComplexFloat> c2(wid2, input.ElemNC(input.Length()-1)); temp1.ReadFromVector(c2, input.Length()); } } Vector<ComplexFloat> temp2; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp2 = Vector<ComplexFloat>(fftlength); for(int i=0; i<filter.Length(); i++) { temp2[i] = filter[i]; } } else if(b == Circular) { temp2 = Vector<ComplexFloat>(input.Length()); for(int i=0; i<filter.Length(); i++) { temp2[i] = filter[i]; } } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } Vector<ComplexFloat> o = IFFT( FFT(temp1) * FFT(temp2) ); if(fa == Same && b != Circular) { o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1); } else if(fa == Valid) { o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length()); } else if(fa == Full && PowerOfTwo) { o = o.Slice(0, input.Length()+filter.Length()-2); } //else if(fa != Full) //{ // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Invalid FilterArea type!"); //} return o; } Vector<ComplexDouble> Conv(Vector<ComplexDouble>& input, Vector<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo) { // Check if the filter size is smaller than the input size if(filter.Length() > input.Length()) { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Filter length should not be larger than the input signal!"); } int fftlength = input.Length()+filter.Length()-1; if(PowerOfTwo) { fftlength = NextPow2(fftlength); } Vector<ComplexDouble> temp1; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp1 = Vector<ComplexDouble>(fftlength); for(int i=0; i<input.Length(); i++) { temp1[i] = input[i]; } } else if(b == Circular) { temp1 = input; } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } if(b == Symmetric) { int wid1 = filter.Length()/2; int wid2 = filter.Length() - wid1 - 1; if(wid1 != 0) { Vector<ComplexDouble> c1 = Reverse( input.Slice(0, wid1-1) ); temp1.ReadFromVector(c1, temp1.Length()-wid1); } if(wid2 != 0) { Vector<ComplexDouble> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) ); temp1.ReadFromVector(c2, input.Length()); } } else if(b == Replicate) { int wid1 = filter.Length()/2; int wid2 = filter.Length() - wid1 - 1; if(wid1 != 0) { Vector<ComplexDouble> c1(wid1, input.ElemNC(0)); temp1.ReadFromVector(c1, temp1.Length()-wid1); } if(wid2 != 0) { Vector<ComplexDouble> c2(wid2, input.ElemNC(input.Length()-1)); temp1.ReadFromVector(c2, input.Length()); } } Vector<ComplexDouble> temp2; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp2 = Vector<ComplexDouble>(fftlength); for(int i=0; i<filter.Length(); i++) { temp2[i] = filter[i]; } } else if(b == Circular) { temp2 = Vector<ComplexDouble>(input.Length()); for(int i=0; i<filter.Length(); i++) { temp2[i] = filter[i]; } } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } Vector<ComplexDouble> o = IFFT( FFT(temp1) * FFT(temp2) ); if(fa == Same && b != Circular) { o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1); } else if(fa == Valid) { o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length()); } else if(fa == Full && PowerOfTwo) { o = o.Slice(0, input.Length()+filter.Length()-2); } //else if(fa != Full) //{ // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Invalid FilterArea type!"); //} return o; } // ============= // Conv2 // ============= Matrix<float> Conv2(Matrix<float>& input, Matrix<float>& filter, FilterArea fa, Border b, bool PowerOfTwo) { return Real( Conv2(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) ); } Matrix<float> Conv2(Vector<float>& filter1, Vector<float>& filter2, Matrix<float>& input, FilterArea fa, Border b, bool PowerOfTwo) { return Real( Conv2(ToComplexFloat(filter1), ToComplexFloat(filter2), ToComplexFloat(input), fa, b, PowerOfTwo) ); } Matrix<double> Conv2(Matrix<double>& input, Matrix<double>& filter, FilterArea fa, Border b, bool PowerOfTwo) { return Real( Conv2(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) ); } Matrix<double> Conv2(Vector<double>& filter1, Vector<double>& filter2, Matrix<double>& input, FilterArea fa, Border b, bool PowerOfTwo) { return Real( Conv2(ToComplexDouble(filter1), ToComplexDouble(filter2), ToComplexDouble(input), fa, b, PowerOfTwo) ); } Matrix<ComplexFloat> Conv2(Matrix<ComplexFloat>& input, Matrix<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo) { // Check if the filter size is smaller than the input size if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns()) { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Filter dimensions should not be larger than the input signal!"); } // Allocate area for input using the border int fftrows = input.Rows()+filter.Rows()-1; int fftcols = input.Columns()+filter.Columns()-1; if(PowerOfTwo) { fftrows = NextPow2(fftrows); fftcols = NextPow2(fftcols); } Matrix<ComplexFloat> temp1; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp1 = Matrix<ComplexFloat>(fftrows, fftcols,0); for(int j=0; j<input.Columns(); j++) { for(int i=0; i<input.Rows(); i++) { temp1.ElemNC(i,j) = input.ElemNC(i,j); } } } else if(b == Circular) { temp1 = input; } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } if(b == Symmetric) { int wid1 = filter.Columns()/2; int wid2 = filter.Columns() - wid1 - 1; if(wid1 != 0) { Matrix<ComplexFloat> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) ); temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1); } if(wid2 != 0) { Matrix<ComplexFloat> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ); temp1.ReadFromMatrix(c2, 0, input.Columns()); } int hei1 = filter.Rows()/2; int hei2 = filter.Rows() - hei1 - 1; if(hei1 != 0) { Matrix<ComplexFloat> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) ); temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0); } if(hei2 != 0) { Matrix<ComplexFloat> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) ); temp1.ReadFromMatrix(r2, input.Rows(), 0); } if(hei1 != 0 && wid1 != 0) { Matrix<ComplexFloat> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) )); temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1); } if(hei1 != 0 && wid2 != 0) { Matrix<ComplexFloat> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) )); temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns()); } if(hei2 != 0 && wid1 != 0) { Matrix<ComplexFloat> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) )); temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1); } if(hei2 != 0 && wid2 != 0) { Matrix<ComplexFloat> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) )); temp1.ReadFromMatrix(q22, input.Rows(), input.Columns()); } } else if(b == Replicate) { int wid1 = filter.Columns()/2; int wid2 = filter.Columns() - wid1 - 1; if(wid1 != 0) { Matrix<ComplexFloat> c1 = input.Slice(0, input.Rows()-1, 0, 0); temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1); } if(wid2 != 0) { Matrix<ComplexFloat> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1); temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns()); } int hei1 = filter.Rows()/2; int hei2 = filter.Rows() - hei1 - 1; if(hei1 != 0) { Matrix<ComplexFloat> r1 = input.Slice(0, 0, 0, input.Columns()-1); temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0); } if(hei2 != 0) { Matrix<ComplexFloat> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1); temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0); } if(hei1 != 0 && wid1 != 0) { Matrix<ComplexFloat> q11( hei1, wid1, input.ElemNC(0,0) ); temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1); } if(hei1 != 0 && wid2 != 0) { Matrix<ComplexFloat> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) ); temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns()); } if(hei2 != 0 && wid1 != 0) { Matrix<ComplexFloat> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) ); temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1); } if(hei2 != 0 && wid2 != 0) { Matrix<ComplexFloat> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1)); temp1.ReadFromMatrix(q22, input.Rows(), input.Columns()); } } // Allocate area for the filter Matrix<ComplexFloat> temp2; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp2 = Matrix<ComplexFloat>(fftrows, fftcols,0); for(int j=0; j<filter.Columns(); j++) { for(int i=0; i<filter.Rows(); i++) { temp2.ElemNC(i,j) = filter.ElemNC(i,j); } } } else if(b == Circular) { temp2 = Matrix<ComplexFloat>(input.Rows(), input.Columns(),0); for(int j=0; j<filter.Columns(); j++) { for(int i=0; i<filter.Rows(); i++) { temp2.ElemNC(i,j) = filter.ElemNC(i,j); } } } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } // FFT and inverse FFT Matrix<ComplexFloat> o = IFFT2( FFT2(temp1)*FFT2(temp2) ); // Based on the filter area, crop the border. if(fa == Same && b != Circular) { o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1); } else if(fa == Valid) { o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns()); } else if(fa == Full && PowerOfTwo) { o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2); } //else if(fa != Full) //{ // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Invalid FilterArea type!"); //} return o; } Matrix<ComplexFloat> Conv2(Vector<ComplexFloat>& filter1, Vector<ComplexFloat>& filter2, Matrix<ComplexFloat>& input, FilterArea fa, Border b, bool PowerOfTwo) { // Check if the filter sizes is smaller than the input size if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns()) { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns"); } Vector<ComplexFloat> *cols = new (std::nothrow) Vector<ComplexFloat>[input.Columns()]; Utility::CheckPointer(cols); for(int i=0; i<input.Columns(); i++) { cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo); } int rlength = cols[0].Length(); Vector<ComplexFloat> *rows = new (std::nothrow) Vector<ComplexFloat>[rlength]; Utility::CheckPointer(rows); for(int i=0; i<rlength; i++) { rows[i] = Vector<ComplexFloat>(input.Columns()); for(int j=0; j<input.Columns(); j++) { rows[i][j] = cols[j][i]; } } delete [] cols; for(int i=0; i<rlength; i++) { rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo); } int clength = rows[0].Length(); Matrix<ComplexFloat> temp(rlength, clength); for(int i=0; i<rlength; i++) { for(int j=0; j<clength; j++) { temp.ElemNC(i,j) = rows[i][j]; } } delete [] rows; return temp; } Matrix<ComplexDouble> Conv2(Matrix<ComplexDouble>& input, Matrix<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo) { // Check if the filter size is smaller than the input size if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns()) { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Filter dimensions should not be larger than the input signal!"); } // Allocate area for input using the border int fftrows = input.Rows()+filter.Rows()-1; int fftcols = input.Columns()+filter.Columns()-1; if(PowerOfTwo) { fftrows = NextPow2(fftrows); fftcols = NextPow2(fftcols); } Matrix<ComplexDouble> temp1; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp1 = Matrix<ComplexDouble>(fftrows, fftcols,0); for(int j=0; j<input.Columns(); j++) { for(int i=0; i<input.Rows(); i++) { temp1.ElemNC(i,j) = input.ElemNC(i,j); } } } else if(b == Circular) { temp1 = input; } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } if(b == Symmetric) { int wid1 = filter.Columns()/2; int wid2 = filter.Columns() - wid1 - 1; if(wid1 != 0) { Matrix<ComplexDouble> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) ); temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1); } if(wid2 != 0) { Matrix<ComplexDouble> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ); temp1.ReadFromMatrix(c2, 0, input.Columns()); } int hei1 = filter.Rows()/2; int hei2 = filter.Rows() - hei1 - 1; if(hei1 != 0) { Matrix<ComplexDouble> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) ); temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0); } if(hei2 != 0) { Matrix<ComplexDouble> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) ); temp1.ReadFromMatrix(r2, input.Rows(), 0); } if(hei1 != 0 && wid1 != 0) { Matrix<ComplexDouble> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) )); temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1); } if(hei1 != 0 && wid2 != 0) { Matrix<ComplexDouble> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) )); temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns()); } if(hei2 != 0 && wid1 != 0) { Matrix<ComplexDouble> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) )); temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1); } if(hei2 != 0 && wid2 != 0) { Matrix<ComplexDouble> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) )); temp1.ReadFromMatrix(q22, input.Rows(), input.Columns()); } } else if(b == Replicate) { int wid1 = filter.Columns()/2; int wid2 = filter.Columns() - wid1 - 1; if(wid1 != 0) { Matrix<ComplexDouble> c1 = input.Slice(0, input.Rows()-1, 0, 0); temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1); } if(wid2 != 0) { Matrix<ComplexDouble> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1); temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns()); } int hei1 = filter.Rows()/2; int hei2 = filter.Rows() - hei1 - 1; if(hei1 != 0) { Matrix<ComplexDouble> r1 = input.Slice(0, 0, 0, input.Columns()-1); temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0); } if(hei2 != 0) { Matrix<ComplexDouble> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1); temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0); } if(hei1 != 0 && wid1 != 0) { Matrix<ComplexDouble> q11( hei1, wid1, input.ElemNC(0,0) ); temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1); } if(hei1 != 0 && wid2 != 0) { Matrix<ComplexDouble> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) ); temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns()); } if(hei2 != 0 && wid1 != 0) { Matrix<ComplexDouble> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) ); temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1); } if(hei2 != 0 && wid2 != 0) { Matrix<ComplexDouble> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1)); temp1.ReadFromMatrix(q22, input.Rows(), input.Columns()); } } // Allocate area for the filter Matrix<ComplexDouble> temp2; if(b == ZeroPad || b == Symmetric || b == Replicate) { temp2 = Matrix<ComplexDouble>(fftrows, fftcols,0); for(int j=0; j<filter.Columns(); j++) { for(int i=0; i<filter.Rows(); i++) { temp2.ElemNC(i,j) = filter.ElemNC(i,j); } } } else if(b == Circular) { temp2 = Matrix<ComplexDouble>(input.Rows(), input.Columns(),0); for(int j=0; j<filter.Columns(); j++) { for(int i=0; i<filter.Rows(); i++) { temp2.ElemNC(i,j) = filter.ElemNC(i,j); } } } else { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Invalid Border type!"); } // FFT and inverse FFT Matrix<ComplexDouble> o = IFFT2( FFT2(temp1) * FFT2(temp2) ); // Based on the filter area, crop the border. if(fa == Same && b != Circular) { o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1); } else if(fa == Valid) { o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns()); } else if(fa == Full && PowerOfTwo) { o = o.Slice(0, input.Rows()+filter.Rows()-2, 0, input.Columns()+filter.Columns()-2); } //else if(fa != Full) //{ // cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; // Utility::RunTimeError("Invalid FilterArea type!"); //} return o; } Matrix<ComplexDouble> Conv2(Vector<ComplexDouble>& filter1, Vector<ComplexDouble>& filter2, Matrix<ComplexDouble>& input, FilterArea fa, Border b, bool PowerOfTwo) { // Check if the filter sizes is smaller than the input size if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns()) { cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns"); } Vector<ComplexDouble> *cols = new (std::nothrow) Vector<ComplexDouble>[input.Columns()]; Utility::CheckPointer(cols); for(int i=0; i<input.Columns(); i++) { cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo); } int rlength = cols[0].Length(); Vector<ComplexDouble> *rows = new (std::nothrow) Vector<ComplexDouble>[rlength]; Utility::CheckPointer(rows); for(int i=0; i<rlength; i++) { rows[i] = Vector<ComplexDouble>(input.Columns()); for(int j=0; j<input.Columns(); j++) { rows[i][j] = cols[j][i]; } } delete [] cols; for(int i=0; i<rlength; i++) { rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo); } int clength = rows[0].Length(); Matrix<ComplexDouble> temp(rlength, clength); for(int i=0; i<rlength; i++) { for(int j=0; j<clength; j++) { temp.ElemNC(i,j) = rows[i][j]; } } delete [] rows; return temp; } };